## Compute HAZUS rebuilding costs

pacman::p_load(tidyverse, fst)

rm(list = ls())

source("code/globals.R")


locationfactors <- read.csv(file.path(INPUT_PUBLIC, "HAZUS", "StateData", "Exported", "MeansAdjFactor.txt")) %>%
  select(CountyFIPS, MeansAdjRes)

############## Cost data from Hazus Technical Manual

table6.3 <- data.frame(
  Class = c(rep("Economy", 4), rep("Average", 4), rep("Custom", 4), rep("Luxury", 4)),
  Floors = c(rep(c("1StoryRes1", "2StoryRes1", "3StoryRes1", "SplitLvlRes1"), 4)),
  basecost = c(
    97.61,
    104.04,
    104.04,
    96.69,
    116.66,
    122.75,
    127.94,
    113.66,
    159.51,
    163.95,
    168.69,
    153.15,
    188.84,
    194.94,
    201.09,
    181.61
  ),
  finishedbasementadj = c(
    26.45,
    15.20,
    15.20,
    15.20,
    32.80,
    21.05,
    16.65,
    21.05,
    53.65,
    30.90,
    22.55,
    30.90,
    59.00,
    34.55,
    25.50,
    34.55
  ),
  unfinishedbasementadj = c(
    9.55,
    6.30,
    6.30,
    6.30,
    11.25,
    7.40,
    5.80,
    7.40,
    21.65,
    12.90,
    9.60,
    12.90,
    22.65,
    13.85,
    10.40,
    13.85
  )
)


table6.4 <- data.frame(
  Class = c(rep("Economy", 3), rep("Average", 3), rep("Custom", 3), rep("Luxury", 3)),
  GarageSize = c(rep(c("1CarGarage", "2CarGarage", "3CarGarage"), 4)),
  avgaddlcost = c(
    18676,
    29263,
    39580,
    19410,
    30224,
    40768,
    21577,
    34319,
    46713,
    25105,
    40136,
    54821
  )
) %>%
  rbind(data.frame( ### See page 6-8 Of Hazus technical manual: Garage cost is assumed to be 0 when it is a carport or there is no garage
    Class = c(rep("Economy", 2), rep("Average", 2), rep("Custom", 2), rep("Luxury", 2)),
    GarageSize = c(rep(c("CarPort", "NoGarage"), 4)),
    avgaddlcost = rep(0, 8)
  ))

table6.5 <- data.frame(
  IncomeRatio = c(
    0,
    0.5,
    0.85,
    1.25,
    2.0
  ),
  Luxury = c(
    0,
    0,
    0,
    0,
    1
  ),
  Custom = c(
    0,
    0,
    0.25,
    1,
    0
  ),
  Average = c(
    0,
    0.25,
    0.75,
    0,
    0
  ),
  Economy = c(
    1,
    0.75,
    0,
    0,
    0
  )
)


a <- read.csv(file.path(INPUT_PUBLIC, "HAZUS", "StateData", "Exported", "BlockInfo.txt")) %>%
  select(
    -OBJECTID, -Tract, -CenLat, -CenLongit,
    -Pct1to2StryRes3, -Pct3to4StryRes3, -Pct5StryplusRes3,
    -PctLowRiseOther, -PctMidRiseOther, -PctHighRiseOther
  ) %>%
  mutate(
    ShareLuxury = case_when(
      IncomeRatio >= 0 & IncomeRatio < 2.0 ~ 0,
      IncomeRatio >= 2.0 ~ 1
    ),
    ShareCustom = case_when(
      IncomeRatio >= 0 & IncomeRatio < 0.85 ~ 0,
      IncomeRatio >= 0.85 & IncomeRatio < 1.25 ~ 0.25,
      IncomeRatio >= 1.25 & IncomeRatio < 2.0 ~ 1,
      IncomeRatio >= 2.0 ~ 0
    ),
    ShareAverage = case_when(
      IncomeRatio >= 0 & IncomeRatio < 0.5 ~ 0,
      IncomeRatio >= 0.5 & IncomeRatio < 0.85 ~ 0.25,
      IncomeRatio >= 0.85 & IncomeRatio < 1.25 ~ 0.75,
      IncomeRatio >= 1.25 ~ 0
    ),
    ShareEconomy = case_when(
      IncomeRatio >= 0 & IncomeRatio < 0.5 ~ 1,
      IncomeRatio >= 0.5 & IncomeRatio < 0.85 ~ 0.75,
      IncomeRatio >= 0.85 ~ 0
    )
  ) %>%
  mutate(
    across(.cols = starts_with("Pct"), .fns = ~ .x / 100) # convert HAZUS-provided percentages to 0-1 instead of 0-100
  ) %>%
  mutate(
    class_shares_check = ShareLuxury + ShareCustom + ShareAverage + ShareEconomy,
    floors_shares_check = Pct1StoryRes1 + Pct2StoryRes1 + Pct3StoryRes1 + PctSplitLvlRes1,
    garage_shares_check = Pct1CarGarage + Pct2CarGarage + Pct3CarGarage + PctCarPort + PctNoGarage
  )

summary(a)
## Make sure shares add to 1
stopifnot(sum(a$class_shares_check != 1) == 0 & sum(a$floors_shares_check != 1) == 0 & sum(a$garage_shares_check != 1) == 0)
## Confirm that all CA Census Blocks are assumed to have 8% of homes with basements, then make a constant to keep track of that value
stopifnot(sum(a$PctWithBasemnt != 0.08) == 0)
pctwithbasement <- 0.08

############################ MAIN STRUCTURE REBUILD COSTS

maincost <- a %>%
  select(
    CensusBlock, IncomeRatio, Pct1StoryRes1, Pct2StoryRes1, Pct3StoryRes1, PctSplitLvlRes1,
    ShareLuxury, ShareCustom, ShareAverage, ShareEconomy
  ) %>%
  pivot_longer(
    cols = starts_with("Share"),
    names_to = "Class",
    names_prefix = "Share",
    values_to = "ClassShare"
  ) %>%
  pivot_longer(
    cols = starts_with("Pct"),
    names_to = "Floors",
    names_prefix = "Pct",
    values_to = "FloorShare"
  ) %>%
  merge(table6.3, by = c("Class", "Floors"), all.x = TRUE) %>%
  mutate(
    WeightedMainCostComponent = ClassShare * FloorShare * basecost,
    WeightedBasementCostComponent = ClassShare * FloorShare * pctwithbasement * finishedbasementadj
  )

m <- maincost %>%
  group_by(CensusBlock) %>%
  summarize(
    IncomeRatio = mean(IncomeRatio),
    obs = n(),
    SumClassShares = sum(ClassShare),
    SumFloorShares = sum(FloorShare),
    SumClassByFloorShares = sum(ClassShare * FloorShare),
    MainCostPerSqFt = sum(WeightedMainCostComponent),
    BasementCostPerSqFt = sum(WeightedBasementCostComponent)
  )

############################ GARAGE REBUILD COSTS

garagecost <- a %>%
  select(
    CensusBlock, IncomeRatio, Pct1CarGarage, Pct2CarGarage, Pct3CarGarage, PctCarPort, PctNoGarage,
    ShareLuxury, ShareCustom, ShareAverage, ShareEconomy
  ) %>%
  pivot_longer(
    cols = starts_with("Share"),
    names_to = "Class",
    names_prefix = "Share",
    values_to = "ClassShare"
  ) %>%
  pivot_longer(
    cols = starts_with("Pct"),
    names_to = "GarageSize",
    names_prefix = "Pct",
    values_to = "GarageSizeShare"
  ) %>%
  merge(table6.4, by = c("Class", "GarageSize"), all.x = TRUE) %>%
  mutate(WeightedGarageCostComponent = ClassShare * GarageSizeShare * avgaddlcost)

g <- garagecost %>%
  group_by(CensusBlock) %>%
  summarize(
    IncomeRatio = mean(IncomeRatio),
    obs = n(),
    SumClassShares = sum(ClassShare),
    SumGarageSizeShares = sum(GarageSizeShare),
    SumClassByFloorShares = sum(ClassShare * GarageSizeShare),
    GarageCostTotal = sum(WeightedGarageCostComponent)
  )

############################ PUT EVERYTHING TOGETHER

stopifnot(nrow(m) == nrow(a) & nrow(g) == nrow(a))

costs <- merge(m[, c("CensusBlock", "IncomeRatio", "MainCostPerSqFt", "BasementCostPerSqFt")],
  g[, c("CensusBlock", "GarageCostTotal")],
  by = "CensusBlock", all.x = TRUE
) %>%
  mutate(BldgCostPerSqFt = MainCostPerSqFt + BasementCostPerSqFt) %>%
  select(-MainCostPerSqFt, -BasementCostPerSqFt) %>%
  mutate(CountyFIPS = substr(CensusBlock, 1, 4)) %>%
  ### Apply R.S. Means Location Adjustment Factors
  merge(locationfactors, by = "CountyFIPS", all.x = TRUE) %>%
  mutate(
    BldgCostPerSqFt = BldgCostPerSqFt * MeansAdjRes,
    GarageCostTotal = GarageCostTotal * MeansAdjRes
  ) %>%
  ### Calculate contents losses and temporary housing/disruption costs
  mutate(
    ContentsPerSqFt = BldgCostPerSqFt * 0.5, ### Contents -- Page 6-9 of the HAZUS Technical Manual Says they use 50% of rebuilding cost
    TempHousingPerSqFtPer24Months = 0.83 * 24 * MeansAdjRes, ### Temporary housing and disruption costs Page 6-14 of HAZUS Technical manual (Assuming 24 months of temp housing). Note I apply R.S. Means Location Factors here too, as a rough way of capturing high rental costs in CA.
    DisruptionPerSqFt = 0.99 * MeansAdjRes
  )

############################# SAVE A BLOCK-LEVEL COST TABLE

costs %>%
  select(
    CensusBlock, BldgCostPerSqFt, GarageCostTotal,
    ContentsPerSqFt, TempHousingPerSqFtPer24Months,
    DisruptionPerSqFt
  ) %>%
  write_fst(file.path(WORKING, "hazuscosts.fst"))







## Save mean and median rebuild costs -----
costs %>%
  ungroup() %>%
  summarize(
    MeanLossesPerSqFt = mean(BldgCostPerSqFt + ContentsPerSqFt + TempHousingPerSqFtPer24Months + DisruptionPerSqFt, na.rm = T),
    MeanGarageTotal = mean(GarageCostTotal, na.rm = T),
    MedianLossesPerSqFt = median(BldgCostPerSqFt + ContentsPerSqFt + TempHousingPerSqFtPer24Months + DisruptionPerSqFt, na.rm = T),
    MedianGarageTotal = median(GarageCostTotal, na.rm = T),
    p25LossesPerSqFt = quantile(BldgCostPerSqFt + ContentsPerSqFt + TempHousingPerSqFtPer24Months + DisruptionPerSqFt, probs = 0.25, na.rm = T),
    p25GarageTotal = quantile(GarageCostTotal, probs = 0.25, na.rm = T),
    p75LossesPerSqFt = quantile(BldgCostPerSqFt + ContentsPerSqFt + TempHousingPerSqFtPer24Months + DisruptionPerSqFt, probs = 0.75, na.rm = T),
    p75GarageTotal = quantile(GarageCostTotal, probs = 0.75, na.rm = T),
    p90LossesPerSqFt = quantile(BldgCostPerSqFt + ContentsPerSqFt + TempHousingPerSqFtPer24Months + DisruptionPerSqFt, probs = 0.90, na.rm = T),
    p90GarageTotal = quantile(GarageCostTotal, probs = 0.90, na.rm = T)
  ) %>%
  write.csv(file.path(WORKING, "loss_info_sample.csv"))
